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1 Introduction 



The collision of particles with the internal degrees of freedom are inelastic in general. 
The inelastic collisions are abundant in naturei . Examples can be seen in collisions 
of atoms, molecules, elastic materials, balls in sports, and so on. The study of 
inelastic collisions will be able to be widely accepted as one of fundamental subjects 
in physics, because they are almost always discussed in textbooks of elementary 
classical mechanics! . 

Recent extensive interest in granular mater ialsi makes physicists to recognize 
fundamental roles of inelastic collisions. In fact, granules consists of macroscopic 
dissipative particles. Therefore, the decision of interaction among particles is obvi- 
ously important. We believe that static interactions among granular particles can 
be described by the theory of elasticityH' i' i i . For example, the normal compres- 
sion may be described by the Hertzian contact forcei and the shear force may be 
represented by the Mindline forcell . The dynamical part related to the dissipation, 
however, cannot be described by any reliable physical theory. Thus, the distinct 
element method0 which is one of the most popular models to simulate collections 
of granular particles contains some dynamical undetermined parameters. In other 
words, to determine such the parameters is important for both granular physics and 
fundamental physics. 

The normal impact of macroscopic materials is characterized by the coefficient 
of restitution (COR) defined by 

e=-Vr/Vi, (1) 
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where Vi and Vr are the relative velocities of incoming and outgoing particles re- 
spectively. COR e had been believed to be a material constant, since the classical 
experiment by Newton0 . In general, however, experiments show that COR for 
three dimensional materials is not a constant even in approximate sense but de- 
pends strongly on the impact velocityS' 111 El . 

The origin of the dissipation in inelastic collisions is the transfer of the kinetic 
energy of the center of mass into the internal degrees of freedom during the impacts. 
Systematic theoretical investigations of the impact have begun with the paper by 
Kuwabara and Kono0 . Taking into account the viscous motion among the internal 
degrees of freedom, they derived the equation of the macroscopic deformation. Later, 
Brilliantov et a/.lll and Morgado and Oppenheimlli derived the identical equation 
to eq.(0). In particular, the derivation by Morgado and Oppenheim0 is based on 
the standard technique of nonequilibrium statistical mechanics to extract the slow 
mode among the fast many modes which can be regarded as the thermal reservoir. 
Furthermore, Brilliantov et al. compared their theoretical results with experimental 
results0 . Thus, the quasi-static theory has been accepted as reasonable one. 

On the other hand, Gerl and Zippelius0 performed the microscopic simulation 
of the two-dimensional collision of an elastic disk with a wall. Their simulation is 
mainly based on the mode expansion of an elastic disk under the force free boundary 
condition. Then, they solve Hamilton's equation determined by the elastic field and 
the repulsive potential to represent the collision of two disks. Their results show that 
COR decreases with the impact velocity, which strongly depends on Poisson's ratio. 
For high velocity of the impact they demonstrate the macroscopic deformation has 
left after the collision is over. Although it is not easy to discuss the impact with the 
very low impact velocity from their method, their analysis may suggest thenossibil- 
ity of a complicated relation between the quasi-static theory of impactB @ and 
their microscopic simulationEl . Thus, we have to clarify the relation between two 
typical approaches. 

In this paper, we will perform the microscopic simulation of the impact of a two 
dimensional elastic disk with a wall. We introduce two methods of simulation; One 
is based on the lattice model (model A) and another is continuum model (model 
B) which is identical to that by Gerl and Zippeliuslli . Through our simulation, 
we will demonstrate that (i) the effect of temperature (the initial internal motion) 
is important, (ii) COR is suddenly dropped by the plastic deformation which is 
enhanced by the initial temperature, and (iii) the continuum model (model B) does 
not recover the results predicted by the quasi-static theories in the low impact 
velocityBlliS. 

The organization of this paper is as follows. In the next section, we will briefly 
review the outline of quasi- static theory^' El . In section 3, we will explain model 
A and model B which is equivalent to the model by Gerl and Zippeliusta of our 
simulation. In section 4, we will show the result of our simulation and discuss the 
validity of quasi-static theory. In section 5, we discuss our results, in particular, 
about the plastic deformation by the impact and its origin. In section 6, we will 
summarize our result. 
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2 Quasi-static theory 



In this section, we briefly explain the outline of quasi-static theory. One purpose 
of this section is to summarize the two-dimensional version of quasi-static theory 
which may not be mentioned in any articles explicitly. 

At first, let us summarize the three dimensional result, in which the equation of 
the macroscopic deformation is given by 

h = -kh^/^ - -fVhh (2) 

in a collision of two spheres, where the macroscopic deformation h is given hj h = 
Ri + R2 — \ri — with the radius Ri{i = 1,2) and the position of the center of the 
mass Fj of i th particle, h and h are respectively dh/dt and d'^h/dt^. k and 7 are 
unimportant constants. The first term of the right hand side in eq.(0) represents 
the Hertzian contact forced' S B i and the second term is the dissipation due to the 
internal motion. 

The simplest derivation of eq. (0) is that by Brilliantov et al^ , though we also 
check it validity by the alternative methods. Taking into account the limitation of 
the length of this paper, we follow the argument by them. 

The static stress tensor in the two-dimensional linear elastic material can be 
represented by 

a^^'hj = 2fi{eij - 5ijeu/2) + KSijeu (3) 

where fi and K are respectively the shear modulus and the compressional modulus, 
and eij is given by 

2[dxj dxj ^' 

with the displacement field Ui. 

The two dimensional Hertzian contact is given by the relation between 

the macroscopic deformation of the center of mass h and the elastic force F^i as 

h ^{In TTT^ ^ - 1}, (5) 



where Y, a and R are the Young modulus, Poisson's ratio and the radius of the 
disk without deformation, respectively. Equation (^) can be derived from the stress 
tensor with the standard treatment of linear elastic theory. Note that h satisfies 
h = R — yo with the position of the center of mass . 

For small dissipation, as in the textbooks! , the dissipative stress tensor due to 
the viscous motion among internal motions is given by 

= 2r],{e,, - 5,,ea/2) + r^^^^.h (6) 

as in the case of viscous fluid, where kij is the time derivative of e^j, ?7j (i = 1, 2) is 
the viscous constant. 

Brilliantov el al^ assumed that the velocity of deformation field is governed by 
the macroscopic deformation, i.e., iii ~ h{dui/dh). Since in the limit of f j we 
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may replace eq.(|^) by Fei ~ —ivYh/ \n{4:R/h)v3 . Thus, with the aid of the assump- 
tion by BriUiantov et alM , (|) and (D, it is easy to derive the two dimensional 
version of quasi-static theory as 

TTYh _ irYh 
" ~ln{4R/h) ~ ln{AR/h)' ^ ' 

where A is an unimportant constant. This result can be derived by various other 
method. Thus, we will compare the result of our simulation with eq.(|^). 

3 Our models 

Let us explain the details of our model. In both models, the wall exists at ?/ = 0, 
and the center of mass keeps the position at x = 0. The disk approaches from y > 
region and is rebounded by the wall. 

3.1 Model A 



Figure 1: A schematic figure of a disk used in model A. 



The disk in model A consists of some mass points (with the mass m) on the 
triangular lattice. All the mass points are combined with linear springs with the 
spring constant k. In the limit of large number of the mass points, this disk corre- 
sponds to the continuum circular disk with the Young's modulus Y = 2k/\/3 and 
Poisson's ratio l/slli. The position of each mass point of model A is governed by 
the following equation: 

= E(rfo - K - r.Dr^ + mVoe'^y" (8) 

oft j^—i I 
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where do is the lattice constant, is the position of the nearest neighbor mass points 
of Tp, m is the mass of the mass points, Up is the y coordinate of r^, and y is the 
unit vector of y direction. Note that the directional projection of the linear spring 
force in model A can cause the nonlinear deformation. The wall potential is given 
by VqC""^ , where Vq = a/2 and a = IOO/c/q for model A. The number of the mass 
points is fixed to 1459 in model A, since the rough evaluation of convergence of the 
results has been checked in this model. 



3.2 Model B 

In this subsection, we introduce model B which is originally proposed by Gerl and 
ZippeliusEl . Although the details of this model can be checked in their paper, 
we present the minimum description of this model to understand the setup of our 
simulation. 

Gerl and Zippeliuslii analyze Hamilton's equation ; 

• _ dH ■ _ OH 

under the Hamiltonian 

H=^ + + l^Mul,Ql,) + Vo f_'" #e-^(^'*). (10) 



n,i 



Here Qn,i is the expansion coefficient of the 2D elastic deformation field in the polar 
coordinate u = (m^, u^) 



{ur{r,(j)),u^{r,(j))) = ^Q„,z«''(r)cosn0,M^''(r)sinn0), (11) 

n.l 



where uy{r)R = An,i h nBn,i — and u^' {r)R = -nAn,i- ^ 

Bn I — ^ with the radius of the disk and the Bessel function of the n— th order 
dr 



Jn{,x). Here k'^ i = kn,iy'2{l + cr)/(l — a^) and kn^i is the solution of 

(1 - Or2)(l - n^)KK'^X,^i{K)Jn-l{K) + K^[k^ - 2n{n + 1)(1 - a)]Jn{K)Jn{K) 

+ (1 - a)[K^ - (1 - - n^)n][Un-i{t^)Jn{f^') + Jn-l{'^')Jni'^)] = (12) 

with Poisson's ratio a, n = knjR and = k'^ iR, which is given by the boundary 
condition. Thus, for fixed n there are infinitely many solutions kn,i and uJn,i = 



kn,iyY/ {p(l — cr^)} numbered by / = 0, 1, ■ ■ ■ , 00. An^i and Bn,i are determined by 

^AdJnik'^lR) Jn{k'nlR), , , 

+ni?„(l - a)[- - -'^\ = (13) 
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and /o^(irr{M"''^ + M0'' } = . Pn^i is the canonical momentum. y{4>,t) is the shape 
of the elastic disk in the polar coordinate; 

y{4>, t) = yo(t) + Qn,i{Cn,i cos(n0) cos - Sn,i sin(n0) sin 0) (14) 

n,l 

with the position of the center of mass yo{t) and constants Cn,i and Sn,i determined 
by the maximal radial and tangential displacement at the edge of the disk as Cn,i = 
u^'\R) and Sn,i = u^\R). M is the mass of the disk, and the momentum of the 
center of the mass po = Myo satisfies po = —{dH/dyo) , Vq and a are parameters to 
express the strength of the wall potential. 

For the simulation of a pair of identical disks, they extrapolate the results of 
their simulation to a — » oo and N ^ oo with the total number of modes A^. We 
only adopt = 1189 (n < 50 and k„ < 50)or N = 437 {n < 30 and k„ < 30), 
Vo = a/2 and a = 500/i? with the radius of the disk R. 

3.3 Parameters in both models 

For the sake of simplicity and comparison between two different models, we only 
simulate the case of Poisson's ratio a = 1/3. The numerical scheme of the integration 
of model A is the classical fourth order Runge-Kutta method with At = 1.6 x 
10~^^Jm/ K. For model B, we adopt the fourth order symplectic integral method 

with At = 5.0 X 10~'^-R/c with c = \Jy/ p for model B where Y is Young's modulus 
and p is the density. In both models, we have checked the conservation of the total 
energy. 

We also investigate the impact with the finite temperature. The temperature is 
introduced as follows: In model A, we prepare the Maxwellian for the initial velocity 
distribution of mass points, where the positions of all mass points are located at their 
equilibrium positions. From the variance of the Maxwellian we can introduce the 
temperature as a parameter. To perform the simulation, we prepare 10 independent 
samples obeying the Maxwellian with the aid of normal random number. In model B, 
we prepare samples in which the absolute value of each mode satisfies equipartition 
law exactly. The sign of each mode is assumed to be at random with the aid of 
the uniform random number. /^From the equipartition law we can introduce the 
temperature parameter of simulation, too. 

The summary of differences between model A and B is as follows: (i) All of 
the mass points in model A interact with the wall but, in model B, only exterior 
boundary has the influence of the potential as in (p!0|). (ii) Model A can have 
nonlinear deformations, but model B is based on the theory of linear elasticity, (iii) 
Model A can express some plastic deformations, but model B cannot, (iv) Model A 
has the six folds symmetry but model B has the rotational symmetry. 

4 Results 

Now, let us explain the details of the result of our simulation. In the first subsection, 
we will introduce the result at T = and in the second subsection, we will show the 
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result at finite T. 



4.1 Simulation at T = 



At first 
condition at T 



we carry out the simulation of model A and model B with the initial 
{i.e. no internal motion). Figure [4.1| is the plot of the COR 







against the impact velocity for both model A and model B. For model B, we show 
the results of 437 modes and 1189 modes which clearly demonstrates the convergence 
of the result for the number of modes. When impact velocity Vi is larger than 0.1c 
with c = \Jy/ p, the value of COR of model A is almost identical to that of model 
B. Each line decreases smoothly as impact velocity increases. At present, we do 
not know the reason why the significant difference between two models exists at low 
impact velocity. 
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Figure 2: Coefficient of restitution for normal collision of the Model A and 
Model B as a function of impact velocity, where c = \/Y/ p with the Young's 
modulus Y and the density p. 



Second, we investigate the force acting on the center of mass of the disk caused 
by the interaction with the wall in model B. In the limit of f , — we expect that the 
Hertzian contact theory can be usedi' B . The small amount of transfer from the 
translational motion to the internal motion is the macroscopic dissipation. Thus, we 
can check the validity of quasi-static approache from our simulation by the 

difference between the observed force acting on the center of mass and the Hertzian 
contact force. 
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Figure 3: The comparison of the Hertzian force in eq.(|5|) with our simulation 
at Vi = 0.01c (a) and Vi = O.OOlc(b) at T = in model B. 



If h is given, we can calculate the elastic force by solving eq.(|^) numerically. 
Figure 2 is the comparison with our simulation in model B (1189 modes) and the 
Hertzian contact theory which is given by the solid lines. The result of our 
simulation at the impact velocity Vi = 0.01c with c = ^JyJp shows the beautiful 
hysteresis as suggested in the simulation at Vi = 0.1c in ref.lll . This means the 
compression and rebound are not symmetric. The hysteresis curve is still self-similar 
even at Vi = 0.04c but the loop becomes noisy at Vi = 0.1c. 

For very low impact velocity Vt = 0.001c, the hysteresis loop disappears but 
the total force observed in our simulation is almost a linear function of h which is 
deviated from the Hertzian contact theory and quasi-static theory (|^. In particular, 
the turning point at F = is deviated from the Hertzian curve (the solid line). 
This deviation is clearly contrast to the quasi static theory, because the dissipative 
force in the theory in eqs.(0) and (|^ must be zero at the turning point at which 
h = should satisfy. This tendency is invariant even for the simulation of model 
A, though the data becomes noisy. The linearity of the total repulsion force is not 
surprising, because e""^*-''^'*^ in the potential term in eq.(|lOD can be expanded by 
series of Qn,i for very slow impact. Although we cannot judge whether the model 
itself is not appropriate for slow impact or the quasi-static theory is wrong, our result 
clearly means that the validity of the quasi-static theory cannot be supported by our 
microscopic simulation. However, the validity of the contact time r in the impact 
evaluated by the quasi-static theory 111 can be evaluated as r ~ {ttR/ c) -y/ln(4c/ Vi) 
has been confirmed by the results of our simulation of model A (Fig. 4). Thus, at 
least, the relation between dynamical impact theory and quasi-static elastic theory 
is not trivial at present. 
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Figure 4: The plot of contact time versus the impact velocity. R represents 
the radius of the disk, in which R = 40,60,70 and 80 correspond to the number 
of mass points 5815, 13057, 17761 and 23233, respectively. The dotted line is 
fitting curve based on the quasi-static theory. 



4.2 Simulation at finite T 

Now, let us show the results of our simulation at finite T which has significant 
differences from those at T = in both low and large impact velocities. In this 
sense, we have much room to study this process at finite T systematically. 

For small impact velocity, COR at finite T becomes larger than that at T = in 
both models. In some trials COR becomes larger than 1. It is an interesting result 
to extract work of this system from thermodynamical point of view. The details of 
the temperature effect at the slow impact will be reported elsewhere. 

For large impact velocity, we do not observe any definite temperature effect in 
model B but we find drastic drop of COR in model A. It seems that COR can be on 
a universal curve when the impact velocity is scaled by the critical velocity above 
which COR drops abruptly (Fig. [4.2| ). The relation between the critical velocity and 
the initial temperature at the intermediate impact velocities is shown in the Fig. |4.2| . 
The critical velocity seems to obey a linear function of T, though the data is not on 
the very slow and the very fast impacts. 



5 discussion 

We investigate what happens in the disk above the critical velocity and find the 
existence of plastic deformation of the disk (Fig.^(a)). Actually, there is no energy 
differences between two configurations in Fig.^(b) which can occur after the strong 
compression during the impact but cannot be released after the impact is over. It 
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Figure 5: The relation the coefficient of restitution and the impact velocity 
rescaled by the critical velocity for each temperature. Curves are plotted in 
the log-log scale. The temperature is scaled by Tq = mc? /ks with the mass 
of the mass points m, and Boltzmann constant /cb- Note that the error bars 
are plotted only in the case T/Tq = 0.03 but is the same order even at other 
T. 
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Figure 6: The plot of the initial temperature and the critical velocity causing 
the plastic deformation. Vcr/c = a{T/To) + b is the fitting curve line from the 
data between T/Tq = 0.02 and 0.05. 
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is well known that the plastic deformation causes the drop of COrI . 

Following the description in ref.6, let us explain the dimensional analysis of the 
two-dimensional plastic deformation. From two-dimensional Hertzian law we 



evaluate h ~ al/R where oq = y iFeiR/ (nY*) with the elastic force Fei, the radius 
without contact R and the effective Young's modulus Y* = — o"^) with Poisson's 
ratio a is the contact length in Hertzian lawi . The work for the compression of the 
disk W is W = {l/2)Mvf ~ Jo' dhF^i ~ Jo^daoao/R, where M and Vi are the 
mass of the disk and the impact velocity, respectively, h* and Oq are respectively 
the maximal compression and and the maximal contact length. Here we neglect the 
logarithmic correction and unimportant numerical factors. Introducing the mean 
contact pressure during dynamical loading pa which satisfies pd ~ Fei/ao, W can be 
evaluated by ~ pd{alY / R. From W ~ Mvf we can express al ~ {MvfR/pdY^'^- 
Let us assume that the impact exceeds the yield pressure for the plastic defor- 
mation. In such the case, the deformation during rebound is frozen. Thus, the 
work in a rebound is W ~ F^h* where F,, is the maximal force during the impact. 
From h* ~ F^/Y* and ~ pdal we evaluate W ~ {pdalY/Y*. Substituting the 
expression of Oq into the expression for W and W we obtain the COR as 

1,? w Y-iMVi-yri' 

— 1/3 

where Vr is the rebound velocity. Thus, we expect the law e ~ f j in the collision 
of a plastic deformed disk. The three dimensional version of evaluation which gives 
e ~ f j ^^'^ agrees well with the experiment! . 

Our finding is, however, something new, because (i) the drop of COR is excited 
by the temperature and (ii) COR decreases more rapidly like e ~ v'^'"^ than that 
for the conventional plastic deformation e ~ "v^^^^ in (15). The mechanism how 



to occur the plastic deformation is not clear at present including the linear law in 
Fig.O. 



6 Conclusion 



We have numerically studied the impact of a two dimensional elastic disk with the 
wall with the aid of model A and model B. The result can be summarized as (i) The 
coefficient of restitution (COR) decreases with the impact velocity, (ii) The result of 
our simulation is not consistent with the result of the two-dimensional quasi-static 
theory. For large impact velocity, there is hysteresis in the deformation of the center 
of mass. For small velocity, there remains the inelastic force even at = 0. (iii) 
There are drastic effects of temperature in both small and large impact velocity, 
(iv) In particular, for large impact velocity of model A, we have found the abrupt 
drop of COR above the critical impact velocity by the plastic deformation. The 
critical velocity of the plastic deformation seems to obey a simple linear function of 
temperature. 

We believe that this preliminary report is meaningful to recognize that physicists 
have poor understanding of such the fundamental process of elementary mechanics. 



11 



(a) 




Figure 7: (a) Plastic deformation of model A with Vi = 0.22 at T = 0.03. 

The solid circle represents the initial circle. The cross points are positions of 
the mass points after the collision, (b) Two configurations are energetically 
equivalent. 



We hope that this letter will invite a lot of interest in the impact from various view 
points. We, at least, have a plan to study three dimensional impacts to clarify the 
relation among the microscopic simulation, experiments and the quasi-static elastic 
theory. 
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